/** 
 * \file gdgemv.cu
 * \author Kyle E. Niemeyer
 * \date 10/04/2011
 *
 * Based on "dgemv.f" from BLAS. Separate functions for dgemv with tranpose of matrix
 * and without tranpose.
 *
 */

////////////////////////////////////////////////////////////////////////

/** gdgemvn performs matrix-vector operations, no transpose. BLAS level two.
 * One of: y = alpha * A * x + beta * y
 *
 * \param[in]   m       number of rows in matrix A
 * \param[in]   n       number of columns in matrix A
 * \param[in]   alpha   scalar alpha
 * \param[in]   A       matrix, dimensions (lda, n)
 * \param[in]   lda     first dimension of A
 * \param[in]   x       vector, dimension (n-1)*|incx|+1
 * \param[in]   incx    increment between elements of x
 * \param[in]   beta    scalar beta
 * \param[in,out] y       array, dimension (n-1)*|incy|+1
 * \param[in]   incy    increment between elements of y
 */
__device__ void gdgemvn ( int m, int n, double alpha, const double *A, int lda, 
                          const double *x, int incx, double beta, double *y, int incy )
{
  
  // check input parameters
  if ( (m <= 0) || (n <= 0) || (lda < max(1,m)) || (incx == 0) || (incy == 0) || ( (alpha == 0.0) && (beta == 1.0) ) ) return;
  
  uint kx;
  if ( incx > 0 ) {
    kx = 1;
  } else {
    kx = 1 - ( n - 1 ) * incx;
  }
  
  uint ky;
  if ( incy > 0 ) {
    ky = 1;
  } else {
    ky = 1 - ( m - 1 ) * incy;
  }
  
  // start operations
  // elements of A are accessed sequentially with one pass through A
  
  // first form y = beta * y
  
  if ( beta != 1.0 ) {
    
    if ( incy == 1 ) {
      
      if ( beta == 0.0 ) {
        for ( uint i = 0; i < m; ++i ) {
          y[i] = 0.0;
        }
      } else {
        for ( uint i = 0; i < m; ++i ) {
          y[i] *= beta;
        }
      }
      
    } else {
      
      uint iy = ky;
      if ( beta == 0.0 ) {
        for ( uint i = 0; i < m; ++i ) {
          y[iy] *= 0.0;
          iy += incy;
        }
      } else {
        for ( uint i = 0; i < m; ++i ) {
          y[iy] *= beta;
          iy += incy;
        }
      }
      
    }
    
  }
  
  if ( alpha == 0.0 ) return;
  
  // form y = alpha * A * x + y
  
  uint jx = kx;
  if ( incy == 1 ) {
    for ( uint j = 0; j < n; ++j ) {
      
      if ( x[jx] != 0.0 ) {
        double temp = alpha * x[jx];
        
        for ( uint i = 0; i < m; ++i ) {
          y[i] += temp * A[i + lda * j];
        }
      }
      
      jx += incx;
    }
  } else {
    for ( uint j = 0; j < n; ++j ) {
      
      if ( x[jx] != 0.0 ) {
        double temp = alpha * x[jx];
        
        uint iy = ky;
        for ( uint i = 0; i < m; ++i ) {
          y[iy] += temp * A[i + lda * j];
          iy += incy;
        }
      }
      
      jx += incx;
    }
  }
  
} // end gdgemvn

////////////////////////////////////////////////////////////////////////

/** gdgemvt performs matrix-vector operations, transpose of matrix. BLAS level two.
 * y = alpha * A' * x + beta * y,
 *
 * \param[in]   m       number of rows in matrix A
 * \param[in]   n       number of columns in matrix A
 * \param[in]   alpha   scalar alpha
 * \param[in]   A       matrix, dimensions (lda, n)
 * \param[in]   lda     first dimension of A
 * \param[in]   x       vector, dimension (n-1)*|incx|+1
 * \param[in]   incx    increment between elements of x
 * \param[in]   beta    scalar beta
 * \param[in,out] y       array, dimension (n-1)*|incy|+1
 * \param[in]   incy    increment between elements of y
 */
__device__ void gdgemvt ( int m, int n, double alpha, const double *A, int lda, 
                          const double *x, int incx, double beta, double *y, int incy )
{
  
  // check input parameters
  if ( (m <= 0) || (n <= 0) || (lda < max(1,m)) || (incx == 0) || (incy == 0) || ( (alpha == 0.0) && (beta == 1.0) ) ) return;
  
  uint kx;
  if ( incx > 0 ) {
    kx = 1;
  } else {
    kx = 1 - ( m - 1 ) * incx;
  }
  
  uint ky;
  if ( incy > 0 ) {
    ky = 1;
  } else {
    ky = 1 - ( n - 1 ) * incy;
  }
  
  // start operations
  // elements of A are accessed sequentially with one pass through A
  
  // first form y = beta * y
  
  if ( beta != 1.0 ) {
    
    if ( incy == 1 ) {
      
      if ( beta == 0.0 ) {
        for ( uint i = 0; i < n; ++i ) {
          y[i] = 0.0;
        }
      } else {
        for ( uint i = 0; i < n; ++i ) {
          y[i] *= beta;
        }
      }
      
    } else {
      
      uint iy = ky;
      if ( beta == 0.0 ) {
        for ( uint i = 0; i < n; ++i ) {
          y[iy] *= 0.0;
          iy += incy;
        }
      } else {
        for ( uint i = 0; i < n; ++i ) {
          y[iy] *= beta;
          iy += incy;
        }
      }
      
    }
    
  }
  
  if ( alpha == 0.0 ) return;
  
  // form y = alpha * A' * x + y
  
  uint jy = ky;
  if ( incx == 1 ) {
    for ( uint j = 0; j < n; ++j ) {
      
      double temp = 0.0;
      
      for ( uint i = 0; i < m; ++i ) {
        temp += A[i + lda * j] * x[i];
      }
      
      y[jy] += alpha * temp;
      jy += incy;
      
    }
  } else {
    for ( uint j = 0; j < n; ++j ) {
      
      double temp = 0.0;
      uint ix = kx;
      
      for ( uint i = 0; i < m; ++i ) {
        temp += A[i + lda * j] * x[ix];
        ix += incx;
      }
      
      y[jy] += alpha * temp;
      jy += incy;
      
    }
  }
  
} // end gdgemvt